We start by loading the packages and functions needed

library(limma)
library(samr)
library(psych)
library(gplots)
library(ggpubr)
library(conover.test)
library(gtools)  #for stars.pval

source("PAINTfunctions.R")

1. Input data and some exploratory analysis

Upon exploratory analysis using all the dataset we do not observe any particular structure.

aaPaint<-read.csv("data/PAINTStudiesBiochemDat.csv")

aaPaint$TimePoint<-factor(aaPaint$TimePoint,levels=c("Day3","Day10"),ordered = T)

str(aaPaint)
## 'data.frame':    140 obs. of  31 variables:
##  $ Study            : chr  "PAINT18" "PAINT18" "PAINT18" "PAINT18" ...
##  $ Group            : chr  "Intervention" "Intervention" "Intervention" "Intervention" ...
##  $ LevelIntervention: chr  "High-18" "High-18" "High-18" "High-18" ...
##  $ TimePoint        : Ord.factor w/ 2 levels "Day3"<"Day10": 2 2 2 2 2 2 2 2 2 2 ...
##  $ PatientID        : int  1821 1822 1823 1824 1825 1826 1827 1828 1830 1851 ...
##  $ Phe              : int  87 86 67 96 45 51 65 53 76 52 ...
##  $ Val              : int  231 211 168 217 93 115 161 119 175 118 ...
##  $ Leu              : int  175 171 120 192 75 109 132 81 164 86 ...
##  $ Iso              : int  71 72 43 74 43 52 50 42 61 49 ...
##  $ Lys              : int  306 340 213 370 134 180 207 195 236 178 ...
##  $ Met              : int  34 50 27 35 17 10 19 27 26 34 ...
##  $ Thr              : int  190 296 434 621 66 309 207 194 290 56 ...
##  $ His              : int  90 120 114 127 62 63 71 75 97 70 ...
##  $ Try              : int  19 20 26 12 19 7 11 21 15 28 ...
##  $ Total.EAA        : int  1203 1366 1212 1744 554 896 923 807 1140 671 ...
##  $ Tyr              : int  23 49 85 45 83 13 23 110 24 92 ...
##  $ Gln              : int  499 724 566 743 370 255 319 365 452 342 ...
##  $ Arg              : int  146 253 154 182 30 117 100 54 162 61 ...
##  $ Cys              : int  17 30 11 14 24 7 11 29 17 20 ...
##  $ Gly              : int  321 474 354 399 229 147 272 284 328 280 ...
##  $ Pro              : int  467 471 285 607 175 262 184 200 252 197 ...
##  $ Glu              : int  132 125 61 167 84 96 62 62 151 150 ...
##  $ Asn              : int  21 28 28 38 32 15 23 52 21 47 ...
##  $ Asp              : int  45 47 17 55 30 28 23 21 40 38 ...
##  $ Ala              : int  396 595 418 885 241 214 210 221 305 307 ...
##  $ Ser              : int  223 285 282 250 146 138 149 151 201 164 ...
##  $ Tau              : int  88 385 116 74 155 40 133 242 153 369 ...
##  $ Orn              : int  416 666 297 522 52 289 320 97 396 82 ...
##  $ Cit              : int  11 11 6 6 14 9 12 20 20 11 ...
##  $ Ammonia          : int  27 74 54 9 67 21 79 67 78 79 ...
##  $ ArgIVIntake      : num  452.7 655.7 254.5 95.9 0 ...
#PCA plots (not taking into account IVArg)
do_PCA_Plot_Nov16((aaPaint[,6:30]),groups = aaPaint$LevelIntervention,legendName = "Arg \n Suppl",scale=F)

do_PCA_Plot_Nov16((aaPaint[,6:30]),groups = aaPaint$TimePoint,legendName = "time",scale=F)

do_PCA_Plot_Nov16((aaPaint[,6:30]),groups = paste(aaPaint$LevelIntervention,aaPaint$TimePoint,sep=""),legendName = "",scale=T)

#PCA plot of day 10 coloured by intervention - no structure appreciated
do_PCA_Plot_Nov16((aaPaint[aaPaint$TimePoint=="Day10",6:30]),groups = aaPaint[aaPaint$TimePoint=="Day10","LevelIntervention"],legendName = "",scale=F)

2. Study correlation of plasma Arginine with IV Arginine intake

Using only the data available for patients on the intervention for PAINT18. We can see from the analysis shown below that while at day 3 there is no correlation between intake and levels of arginine, by day 10 there is a strong correlation (Spearman correlation 0.73 with p value <0.001). It is interesting to see that the patients that were on enteral feeds (zero levels at Day 10 for IV Arginine), show relatively low levels of Arginine. Note patient 1824 died the day after D10 sample.

ArgCorre<-aaPaint[!is.na(aaPaint$ArgIVIntake),]
ArgCorre$TimePoint<-factor(ArgCorre$TimePoint,levels=c("Day3","Day10"),ordered = T)

ggplot(data=ArgCorre,aes(x=Arg,y=ArgIVIntake))+
  geom_point()+ggrepel::geom_label_repel(aes(label=PatientID),box.padding=0.35,point.padding = 0.5,segment.color="grey50")+facet_wrap(~TimePoint,scales="free_x")+
  geom_smooth(aes(x=Arg,y=ArgIVIntake),method='lm')+
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 10))+
  xlab("Arginine levels plasma")+ylab("Arginine IV Intake")+
  stat_cor(aes(x=Arg,y=ArgIVIntake),method="spearman", col="blue",label.y = 775)+
  theme_bw(base_size = 16)

3. Statistical analysis of total essential aminoacids, plasma arginine and amonia

We start by producing some visualisations

aaPaintM<-melt(aaPaint[,-ncol(aaPaint)],id.vars=c("Study","Group","LevelIntervention","TimePoint","PatientID"))
aaPaintM$variable<-as.character(aaPaintM$variable)

aaPaintM_EAA_Arg_Am<-aaPaintM[aaPaintM$variable%in%c("Total.EAA","Arg","Ammonia"),]

ggplot(data=aaPaintM_EAA_Arg_Am)+
  geom_boxplot(aes(x=LevelIntervention, y=value,fill=TimePoint),col="black")+
  geom_point(aes(x=LevelIntervention, y=value,fill=TimePoint,shape=Study),position=position_jitterdodge(dodge.width=0.9),col="black",alpha=0.4,pch=21)+
  facet_wrap(~variable,scales="free")+
  ylab("")+
  scale_fill_manual(values = c("#DCE319FF","#33638DFF"))+
  scale_colour_manual(values = c("#DCE319FF","#33638DFF"))+
  scale_shape_manual(values=c(21,22,23))+
  theme_bw(base_size = 16)+theme(legend.position = "bottom")

aaPaintM_sel<-aaPaintM[aaPaintM$variable%in%c("Total.EAA","Arg","Ammonia","Cit","Orn"),]

ggplot(data=aaPaintM_sel)+
  geom_boxplot(aes(x=LevelIntervention, y=value,fill=TimePoint),col="black")+
  geom_point(aes(x=LevelIntervention, y=value,fill=TimePoint,shape=Study),position=position_jitterdodge(dodge.width=0.9),col="black",alpha=0.4,pch=21)+
  facet_wrap(~variable,scales="free")+
  ylab("")+
  scale_fill_manual(values = c("#DCE319FF","#33638DFF"))+
  scale_colour_manual(values = c("#DCE319FF","#33638DFF"))+
  scale_shape_manual(values=c(21,22,23))+
  theme_bw(base_size = 16)+theme(legend.position = "bottom")

Comparisons for Day 10

Upon data exploration the statistical test chosen to look at the differences at day 10 for each metabolite is a Kruskal-Wallis, corrected for FDR with Benjamini-Hochberg followed by a Conover-Iman test as post-Hoc. Significant aminoacids shown in boxplots with standard significance values.

qqnorm(aaPaint$Arg)#no

shapiro.test(aaPaint$Arg)
## 
##  Shapiro-Wilk normality test
## 
## data:  aaPaint$Arg
## W = 0.7976, p-value = 1.255e-12
qqnorm(aaPaint$Ammonia)#yes

shapiro.test(aaPaint$Ammonia)
## 
##  Shapiro-Wilk normality test
## 
## data:  aaPaint$Ammonia
## W = 0.98934, p-value = 0.362
qqnorm(aaPaint$Total.EAA)#No

shapiro.test(aaPaint$Total.EAA)
## 
##  Shapiro-Wilk normality test
## 
## data:  aaPaint$Total.EAA
## W = 0.9713, p-value = 0.004793
#with the above prelim tests it is better to undertake a non-parametric test
#we'll start by doing a Kruskall Wallis at all the day 10 measurements followed by post-hoc with conover test


aaPaint_D10<-aaPaint[aaPaint$TimePoint=="Day10",]
KW_all_pvals<-apply(aaPaint_D10[,6:30],2,function(x){kruskal.test(x,aaPaint_D10$LevelIntervention)$p.value})

#we correct for FDR and we select the significant for the posthoc

KW_all_pvals_adj<-p.adjust(KW_all_pvals,method = "BH")



namesSig<-names(KW_all_pvals_adj[KW_all_pvals_adj<0.05])#we get 8 significant 
forConover<-aaPaint_D10[,colnames(aaPaint_D10)%in%namesSig]



KW_conoverIman_postHoc<-do.call(rbind,lapply(lapply(forConover,conover.test, g=aaPaint_D10$LevelIntervention,list=T,method="BH",label=T,table=F),as.data.frame))
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 18.9695, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -2.194453 (0.0159)*
## High-18 - Standard_6  : -4.970219 (0.0000)*
## Moderate - Standard_6 : -3.134436 (0.0019)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 10.32, df = 2, p-value = 0.01
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -0.649611 (0.2591)
## High-18 - Standard_6  : -3.092535 (0.0044)*
## Moderate - Standard_6 : -2.747346 (0.0058)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 7.3808, df = 2, p-value = 0.02
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    :  0.064360 (0.4744)
## High-18 - Standard_6  : -2.230323 (0.0219)*
## Moderate - Standard_6 : -2.574716 (0.0185)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 11.8363, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -3.575898 (0.0010)*
## High-18 - Standard_6  : -1.283648 (0.1019)
## Moderate - Standard_6 :  2.541188 (0.0101)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 8.5445, df = 2, p-value = 0.01
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -3.059427 (0.0048)*
## High-18 - Standard_6  : -2.135894 (0.0273)
## Moderate - Standard_6 :  1.009632 (0.1582)
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 20.9383, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    :  1.772715 (0.0405)
## High-18 - Standard_6  :  5.174610 (0.0000)*
## Moderate - Standard_6 :  3.833430 (0.0002)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 19.2275, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -3.725538 (0.0003)*
## High-18 - Standard_6  :  0.568754 (0.2857)
## Moderate - Standard_6 :  4.786731 (0.0000)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 12.9503, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -3.933887 (0.0003)*
## High-18 - Standard_6  : -2.108928 (0.0291)*
## Moderate - Standard_6 :  2.013617 (0.0241)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 7.5947, df = 2, p-value = 0.02
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -1.124172 (0.1325)
## High-18 - Standard_6  : -2.798143 (0.0101)*
## Moderate - Standard_6 : -1.888525 (0.0476)
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 19.5445, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -2.005475 (0.0246)*
## High-18 - Standard_6  :  2.696496 (0.0067)*
## Moderate - Standard_6 :  5.198173 (0.0000)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 15.5167, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -0.364130 (0.3585)
## High-18 - Standard_6  : -3.747581 (0.0003)*
## Moderate - Standard_6 : -3.800377 (0.0005)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 9.2122, df = 2, p-value = 0.01
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -3.153203 (0.0037)*
## High-18 - Standard_6  : -1.373626 (0.0871)
## Moderate - Standard_6 :  1.969532 (0.0399)
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 20.0594, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -0.935407 (0.1765)
## High-18 - Standard_6  :  3.623966 (0.0004)*
## Moderate - Standard_6 :  5.067460 (0.0000)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 16.4153, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -0.034309 (0.4864)
## High-18 - Standard_6  :  3.653762 (0.0004)*
## Moderate - Standard_6 :  4.138753 (0.0002)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 7.7156, df = 2, p-value = 0.02
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate    : -2.759684 (0.0113)*
## High-18 - Standard_6  : -2.381707 (0.0152)*
## Moderate - Standard_6 :  0.441281 (0.3303)
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2
KW_conoverIman_postHoc$AminoAcid<-strsplit2(rownames(KW_conoverIman_postHoc),"[.]")[,1]

This is the table with the post-hoc analysis for the significant aminoacids.

KW_conoverIman_postHoc[KW_conoverIman_postHoc$AminoAcid%in%c("Total"),"AminoAcid"]<-"Total.EAA"

KW_conoverIman_postHoc[,c(6,5,2,4)]
##             AminoAcid           comparisons           T   P.adjusted
## Phe.1             Phe    High-18 - Moderate -2.19445336 1.589113e-02
## Phe.2             Phe  High-18 - Standard_6 -4.97021968 7.721494e-06
## Phe.3             Phe Moderate - Standard_6 -3.13443624 1.937981e-03
## Val.1             Val    High-18 - Moderate -0.64961127 2.591162e-01
## Val.2             Val  High-18 - Standard_6 -3.09253593 4.384765e-03
## Val.3             Val Moderate - Standard_6 -2.74734601 5.824426e-03
## Leu.1             Leu    High-18 - Moderate  0.06436025 4.744405e-01
## Leu.2             Leu  High-18 - Standard_6 -2.23032310 2.188851e-02
## Leu.3             Leu Moderate - Standard_6 -2.57471604 1.847714e-02
## Thr.1             Thr    High-18 - Moderate -3.57589808 9.972381e-04
## Thr.2             Thr  High-18 - Standard_6 -1.28364837 1.019112e-01
## Thr.3             Thr Moderate - Standard_6  2.54118872 1.008235e-02
## Total.EAA.1 Total.EAA    High-18 - Moderate -3.05942756 4.830209e-03
## Total.EAA.2 Total.EAA  High-18 - Standard_6 -2.13589452 2.734474e-02
## Total.EAA.3 Total.EAA Moderate - Standard_6  1.00963286 1.582075e-01
## Arg.1             Arg    High-18 - Moderate  1.77271522 4.048072e-02
## Arg.2             Arg  High-18 - Standard_6  5.17461063 3.578830e-06
## Arg.3             Arg Moderate - Standard_6  3.83343065 2.157958e-04
## Cys.1             Cys    High-18 - Moderate -3.72553847 3.076751e-04
## Cys.2             Cys  High-18 - Standard_6  0.56875407 2.857415e-01
## Cys.3             Cys Moderate - Standard_6  4.78673184 1.524334e-05
## Gly.1             Gly    High-18 - Moderate -3.93388718 3.087180e-04
## Gly.2             Gly  High-18 - Standard_6 -2.10892838 2.910623e-02
## Gly.3             Gly Moderate - Standard_6  2.01361725 2.409668e-02
## Pro.1             Pro    High-18 - Moderate -1.12417239 1.325372e-01
## Pro.2             Pro  High-18 - Standard_6 -2.79814379 1.013419e-02
## Pro.3             Pro Moderate - Standard_6 -1.88852595 4.756429e-02
## Asn.1             Asn    High-18 - Moderate -2.00547530 2.460601e-02
## Asn.2             Asn  High-18 - Standard_6  2.69649638 6.733013e-03
## Asn.3             Asn Moderate - Standard_6  5.19817339 3.474889e-06
## Asp.1             Asp    High-18 - Moderate -0.36413094 3.584710e-01
## Asp.2             Asp  High-18 - Standard_6 -3.74758181 2.862929e-04
## Asp.3             Asp Moderate - Standard_6 -3.80037745 4.814110e-04
## Ser.1             Ser    High-18 - Moderate -3.15320345 3.666462e-03
## Ser.2             Ser  High-18 - Standard_6 -1.37362619 8.713847e-02
## Ser.3             Ser Moderate - Standard_6  1.96953297 3.986740e-02
## Tau.1             Tau    High-18 - Moderate -0.93540725 1.765473e-01
## Tau.2             Tau  High-18 - Standard_6  3.62396628 4.315131e-04
## Tau.3             Tau Moderate - Standard_6  5.06746037 5.512059e-06
## Orn.1             Orn    High-18 - Moderate -0.03430907 4.863679e-01
## Orn.2             Orn  High-18 - Standard_6  3.65376272 3.883769e-04
## Orn.3             Orn Moderate - Standard_6  4.13875353 1.537773e-04
## Cit.1             Cit    High-18 - Moderate -2.75968448 1.134998e-02
## Cit.2             Cit  High-18 - Standard_6 -2.38170764 1.519900e-02
## Cit.3             Cit Moderate - Standard_6  0.44128108 3.302602e-01
write.csv(KW_conoverIman_postHoc,"D10_KW_ConoverImanPH.csv")

Now we generate a plot with the results of the significant aminoacids.

aaPaintM_D10_sig<-aaPaintM[aaPaintM$TimePoint=="Day10" & aaPaintM$variable%in%namesSig,]
colnames(aaPaintM_D10_sig)[6]<-"aminoacid"

KW_conoverIman_postHoc$group1<-strsplit2(KW_conoverIman_postHoc$comparisons," - ")[,1]
KW_conoverIman_postHoc$group2<-strsplit2(KW_conoverIman_postHoc$comparisons," - ")[,2]
ypos<-apply(aaPaint[,namesSig],2,max,na.rm=T)*1.12
KW_conoverIman_postHoc$y.position<-ypos[match(KW_conoverIman_postHoc$AminoAcid,names(ypos))]
KW_conoverIman_postHoc$StatSig<-stars.pval(KW_conoverIman_postHoc$P.adjusted)
KW_conoverIman_postHoc$StatSig[KW_conoverIman_postHoc$StatSig==" "]<-"NS"

KW_conoverIman_postHoc[,c(6,5,2,4)]
##             AminoAcid           comparisons           T   P.adjusted
## Phe.1             Phe    High-18 - Moderate -2.19445336 1.589113e-02
## Phe.2             Phe  High-18 - Standard_6 -4.97021968 7.721494e-06
## Phe.3             Phe Moderate - Standard_6 -3.13443624 1.937981e-03
## Val.1             Val    High-18 - Moderate -0.64961127 2.591162e-01
## Val.2             Val  High-18 - Standard_6 -3.09253593 4.384765e-03
## Val.3             Val Moderate - Standard_6 -2.74734601 5.824426e-03
## Leu.1             Leu    High-18 - Moderate  0.06436025 4.744405e-01
## Leu.2             Leu  High-18 - Standard_6 -2.23032310 2.188851e-02
## Leu.3             Leu Moderate - Standard_6 -2.57471604 1.847714e-02
## Thr.1             Thr    High-18 - Moderate -3.57589808 9.972381e-04
## Thr.2             Thr  High-18 - Standard_6 -1.28364837 1.019112e-01
## Thr.3             Thr Moderate - Standard_6  2.54118872 1.008235e-02
## Total.EAA.1 Total.EAA    High-18 - Moderate -3.05942756 4.830209e-03
## Total.EAA.2 Total.EAA  High-18 - Standard_6 -2.13589452 2.734474e-02
## Total.EAA.3 Total.EAA Moderate - Standard_6  1.00963286 1.582075e-01
## Arg.1             Arg    High-18 - Moderate  1.77271522 4.048072e-02
## Arg.2             Arg  High-18 - Standard_6  5.17461063 3.578830e-06
## Arg.3             Arg Moderate - Standard_6  3.83343065 2.157958e-04
## Cys.1             Cys    High-18 - Moderate -3.72553847 3.076751e-04
## Cys.2             Cys  High-18 - Standard_6  0.56875407 2.857415e-01
## Cys.3             Cys Moderate - Standard_6  4.78673184 1.524334e-05
## Gly.1             Gly    High-18 - Moderate -3.93388718 3.087180e-04
## Gly.2             Gly  High-18 - Standard_6 -2.10892838 2.910623e-02
## Gly.3             Gly Moderate - Standard_6  2.01361725 2.409668e-02
## Pro.1             Pro    High-18 - Moderate -1.12417239 1.325372e-01
## Pro.2             Pro  High-18 - Standard_6 -2.79814379 1.013419e-02
## Pro.3             Pro Moderate - Standard_6 -1.88852595 4.756429e-02
## Asn.1             Asn    High-18 - Moderate -2.00547530 2.460601e-02
## Asn.2             Asn  High-18 - Standard_6  2.69649638 6.733013e-03
## Asn.3             Asn Moderate - Standard_6  5.19817339 3.474889e-06
## Asp.1             Asp    High-18 - Moderate -0.36413094 3.584710e-01
## Asp.2             Asp  High-18 - Standard_6 -3.74758181 2.862929e-04
## Asp.3             Asp Moderate - Standard_6 -3.80037745 4.814110e-04
## Ser.1             Ser    High-18 - Moderate -3.15320345 3.666462e-03
## Ser.2             Ser  High-18 - Standard_6 -1.37362619 8.713847e-02
## Ser.3             Ser Moderate - Standard_6  1.96953297 3.986740e-02
## Tau.1             Tau    High-18 - Moderate -0.93540725 1.765473e-01
## Tau.2             Tau  High-18 - Standard_6  3.62396628 4.315131e-04
## Tau.3             Tau Moderate - Standard_6  5.06746037 5.512059e-06
## Orn.1             Orn    High-18 - Moderate -0.03430907 4.863679e-01
## Orn.2             Orn  High-18 - Standard_6  3.65376272 3.883769e-04
## Orn.3             Orn Moderate - Standard_6  4.13875353 1.537773e-04
## Cit.1             Cit    High-18 - Moderate -2.75968448 1.134998e-02
## Cit.2             Cit  High-18 - Standard_6 -2.38170764 1.519900e-02
## Cit.3             Cit Moderate - Standard_6  0.44128108 3.302602e-01
write.csv(KW_conoverIman_postHoc[,c(6,5,2,4)],"ConoverImanResults.csv")

for(i in namesSig){
paa<-ggplot(data=aaPaintM_D10_sig[aaPaintM_D10_sig$aminoacid==i,],aes(x=LevelIntervention))+
  geom_boxplot(aes(x=LevelIntervention, y=value),col="black")+
  geom_jitter(aes(x=LevelIntervention, y=value),col="black",alpha=0.6,pch=16)+
  facet_wrap(~aminoacid,scales="free")+
  ylab("")+xlab("")+
  theme_classic(base_size = 16)+theme(legend.position = "bottom")



paa2<-paa+stat_pvalue_manual(data=KW_conoverIman_postHoc[KW_conoverIman_postHoc$AminoAcid==i,],label="StatSig",step.increase = c(0,0.1,0.1))
 

ggsave(plot = paa2,filename = paste("FigBiochem/",i,"_Day10_boxplot.png",sep = ""))

print(paa2)

}

Comparisons for day 3

While the first strategy included looking at comparisons between the three groups, as clarified by the clinical team not all the studies are equivalent as some started the intervention after the sample was taken (PAINT 18). Hence Day 3 comparisons will be done in a 3-side non parametric test

Day 3 PAINT 18 data

#for PAINT18
aaPaint_D3_PAINT18<-aaPaint[aaPaint$TimePoint=="Day3" & aaPaint$Study=="PAINT18",]
D3PAINT18_stats_pvals<-sapply(6:30,function(x){
  wilcox.test(x = aaPaint_D3_PAINT18[aaPaint_D3_PAINT18$LevelIntervention=="Standard_6",x],
              y = aaPaint_D3_PAINT18[aaPaint_D3_PAINT18$LevelIntervention=="High-18",x])$p.value})

D3PAINT18_stats<-data.frame("aminoacid"=colnames(aaPaint_D3_PAINT18)[6:30],"pvalue"=D3PAINT18_stats_pvals,"Padj"=p.adjust(D3PAINT18_stats_pvals,method = "BH"))

write.csv(D3PAINT18_stats,"D3PAINT18_CvsIComparison_SuppTable3.csv")

Day 3 PAINT data

aaPaint_D3_PAINT<-aaPaint[aaPaint$TimePoint=="Day3" & aaPaint$Study%in%c("PAINT","PAINT-NH3"),]

D3PAINT_stats_pvals<-sapply(6:30,function(x){
  wilcox.test(x = aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Standard_6",x],
              y = aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Moderate",x])$p.value})

D3PAINT_stats<-data.frame("aminoacid"=colnames(aaPaint_D3_PAINT)[6:30],"pvalue"=D3PAINT_stats_pvals,"Padj"=p.adjust(D3PAINT_stats_pvals,method = "BH"))

D3PAINT_stats$medianControl<-apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Standard_6",6:30],2,median,na.rm=T)
D3PAINT_stats$medianInterv<-apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Moderate",6:30],2,median,na.rm=T)
D3PAINT_stats$IQR_Control<-apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Standard_6",6:30],2,IQR,na.rm=T)
D3PAINT_stats$IQR_Interv<-apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Moderate",6:30],2,IQR,na.rm=T)

write.csv(D3PAINT_stats,"D3PAINT_CvsIComparison_SuppTable2.csv")


apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Standard_6",6:30],2,quantile,na.rm=T)
##         Phe    Val    Leu   Iso    Lys   Met    Thr   His  Try Total.EAA    Tyr
## 0%    59.00 113.00  87.00 13.00 135.00 13.00 135.00  54.0  9.0    788.00  11.00
## 25%   66.25 160.25 131.25 46.00 172.75 22.75 260.25  76.5 14.0   1030.50  32.25
## 50%   88.50 180.00 145.00 50.50 212.50 27.50 325.00  85.5 24.5   1220.00  79.00
## 75%   98.00 204.25 157.00 64.25 275.25 40.25 511.00 106.0 32.5   1427.75 143.25
## 100% 124.00 272.00 204.00 80.00 424.00 50.00 626.00 141.0 41.0   1765.00 657.00
##         Gln   Arg   Cys   Gly    Pro    Glu   Asn    Asp   Ala    Ser   Tau
## 0%   157.00  5.00  5.00 237.0 215.00  60.00 15.00 135.00 197.0  90.00  17.0
## 25%  358.75 11.75 10.75 275.5 321.50 101.25 19.75 172.75 275.5 166.50  42.5
## 50%  442.50 19.50 16.50 351.0 339.00 121.00 24.00 212.50 389.5 225.00  57.0
## 75%  527.00 36.25 22.75 447.5 440.25 177.75 33.25 275.25 428.5 280.75  86.0
## 100% 664.00 61.00 42.00 748.0 492.00 267.00 47.00 424.00 556.0 338.00 187.0
##         Orn   Cit Ammonia
## 0%    22.00 11.00   28.00
## 25%   68.25 12.75   43.50
## 50%   91.00 14.50   52.50
## 75%  120.00 19.50   63.25
## 100% 184.00 27.00   78.00
apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Moderate",6:30],2,quantile,na.rm=T)
##      Phe    Val    Leu    Iso    Lys   Met    Thr   His  Try Total.EAA   Tyr
## 0%    63 109.00  41.00  30.00 136.00 13.00 113.00  55.0 11.0     696.0  17.0
## 25%   74 137.75 115.75  43.75 173.50 22.75 208.75  68.0 16.5     929.0  40.5
## 50%   85 179.50 145.50  51.50 206.00 28.00 267.50  90.0 24.0    1072.0  74.0
## 75%  104 201.25 155.50  58.00 250.75 35.00 293.75 108.5 32.5    1290.5 105.5
## 100% 129 331.00 289.00 103.00 365.00 57.00 463.00 137.0 56.0    1823.0 332.0
##         Gln    Arg  Cys   Gly   Pro   Glu  Asn   Asp   Ala    Ser    Tau   Orn
## 0%   214.00  13.00  3.0 209.0 181.0  74.0 14.0 19.00 183.0 116.00  13.00  60.0
## 25%  348.75  28.50  7.0 274.5 243.5 107.5 21.5 25.75 276.5 153.50  40.25 112.5
## 50%  467.50  40.00 14.5 327.5 276.0 137.5 27.0 30.50 320.0 213.50  77.50 165.5
## 75%  540.25  57.25 23.5 411.0 344.5 171.5 33.0 39.25 459.0 269.75 120.50 227.0
## 100% 676.00 164.00 49.0 510.0 613.0 273.0 78.0 72.00 656.0 390.00 264.00 468.0
##        Cit Ammonia
## 0%    8.00    9.00
## 25%  13.00   38.25
## 50%  14.00   41.00
## 75%  17.25   46.00
## 100% 26.00   64.00

Day 30 PAINT 18 data

#for PAINT18
aaPaintWithD30<-read.csv("data/PAINTStudiesBiochemDat_withD30.csv")
aaPaint_D30_PAINT18<-aaPaintWithD30[aaPaintWithD30$TimePoint=="Day30" & aaPaintWithD30$Study=="PAINT18",]
D30PAINT18_stats_pvals<-sapply(6:30,function(x){
  wilcox.test(x = aaPaint_D30_PAINT18[aaPaint_D30_PAINT18$LevelIntervention=="Standard_6",x],
              y = aaPaint_D30_PAINT18[aaPaint_D30_PAINT18$LevelIntervention=="High-18",x],na.rm=T)$p.value})

D30PAINT18_stats<-data.frame("aminoacid"=colnames(aaPaint_D30_PAINT18)[6:30],"pvalue"=D30PAINT18_stats_pvals,"Padj"=p.adjust(D30PAINT18_stats_pvals,method = "BH"))

write.csv(D30PAINT18_stats,"D30PAINT18_CvsIComparison_SuppTable4.csv")

Below find preparation of some of the figures asked by the medical team for the manuscript:

#boxplots comparing key variables at day 3

aaPaintM<-melt(aaPaint[,-ncol(aaPaint)],id.vars=c("Study","Group","LevelIntervention","TimePoint","PatientID"))
aaPaintM$variable<-as.character(aaPaintM$variable)

aaPaintM_EAA_Arg_Am_Day3<-aaPaintM[aaPaintM$variable%in%c("Total.EAA","Arg","Ammonia","Gln","Cit","Orn") & aaPaintM$TimePoint=="Day3",]


boxplotD3<-ggplot(data=aaPaintM_EAA_Arg_Am_Day3)+
  geom_boxplot(aes(x=LevelIntervention, y=value),col="black",outlier.color = "white")+
  geom_point(aes(x=LevelIntervention, y=value,shape=Study),alpha=0.6)+
  facet_wrap(~variable,scales="free")+
  ylab("Abundance")+ xlab("Intervention")+
  theme_classic(base_size = 12)+theme(legend.position = "bottom")
boxplotD3

#Day 10
aaPaintM_EAA_Arg_Am_Day10<-aaPaintM[aaPaintM$variable%in%c("Total.EAA","Arg","Ammonia","Gln","Cit","Orn") & aaPaintM$TimePoint=="Day10",]

aaPaintM_EAA_Arg_Am_Day10$variable<-factor(aaPaintM_EAA_Arg_Am_Day10$variable,ordered=T,levels=c("Arg","Total.EAA","Ammonia","Orn","Cit","Gln"))
aaPaintM_EAA_Arg_Am_Day10$Study[aaPaintM_EAA_Arg_Am_Day10$Study%in%c("PAINT","PAINT-NH3")]<-"PAINT"

boxplotD10<-ggplot(data=aaPaintM_EAA_Arg_Am_Day10)+
  geom_boxplot(aes(x=LevelIntervention, y=value),col="black",outlier.color = "white")+
  geom_point(aes(x=LevelIntervention, y=value,shape=Study),alpha=0.6)+
  facet_wrap(~variable,scales="free")+
  ylab("Abundance")+ xlab("Intervention")+
  theme_classic(base_size = 12)+theme(legend.position = "bottom")
boxplotD10

ggsave(plot=boxplotD10,filename = "D10_SelectedAminoAcids.png",device = "png",dpi = 300,width=8,units = "in")

ggplot(data=aaPaintM[aaPaintM$variable%in%c("Total.EAA","Arg","Ammonia","Gln","Cit","Orn"),])+
  geom_boxplot(aes(x=LevelIntervention, y=value, fill=TimePoint),col="black",outlier.color = "white")+
  geom_point(aes(x=LevelIntervention, y=value,shape=Study,col=TimePoint,fill=TimePoint),alpha=0.6,position=position_jitterdodge(dodge.width=0.9))+
  scale_fill_manual(values = c("darkgrey","skyblue"))+
  scale_colour_manual(values = c("darkgrey","skyblue"))+
  facet_wrap(~variable,scales="free")+
  ylab("Abundance")+ xlab("Intervention")+
  theme_classic(base_size = 12)+theme(legend.position = "bottom")